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Molecular dynamics simulations of ultrathin water film confined between atomically flat rigid diamond plates 
are described. Films with thickness of one and two molecular diameters are concerned and TIP4P model is 
used for water molecules. Dynamical and equilibrium characteristics of the system for different values of the 
external load and shear force are investigated. An increase of the external load causes the transition of the 
film to a solidlike state. This is manifested in the decreasing of the diffusion constant and in the ordering 
of the liquid molecules into quasidiscrete layers. For two-layer film under high loads molecules also become 
ordered parallel to the surfaces. Time dependencies of the friction force and the changes of its average value 
with the load are obtained. In general, the behaviour of the studied model is consistent with the experimental 
results obtained for simple liquids with spherical molecules. 
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1. Introduction 

Studying of friction on the atomic level is important from fundamental and practical sides. The 
development of nanotechnology and the creation of new devices on its base require a deep under- 
standing of the nanoscale processes leading to wear and friction [HJEJGJH]. I n the nanotribological 
studies there is a phenomenon which has a special place and is called boundary lubrication. It 
corresponds to the situation when rubbing surfaces are separated by an ultrathin (a few molecular 
diameters) lubricant film. Boundary lubrication is important not only for the nanotechnology. It 
is observed under the ordinary conditions, because almost always there is a lubricant between the 
solids (called "the third bodies" by tribologists) - either a specially chosen lubricant film, or wear 
debris produced by sliding, or water adsorbed from air, etc. Boundary lubrication is also important 
in the macromachines at stop/start moments, when the lubricant is squeezed out from the contact 
area and the surfaces come into direct contact [ [TJ [2j |3j [4] . 

As show experiments with surface-force apparatus (SFA) and the computer simulations the 
behaviour of ultrathin liquid film confined between two solid substrates entirely differs from the 
bulk one. The main features of a boundary lubricant are [ [TJ [H [5j [6] : 

1. Liquid molecules confined between two atomically flat surfaces become more ordered and 
tend to organize quasidiscrete layers, where the average local density of the liquid oscillates 
with distance normal to the boundaries. 

2. The mobility of molecules in confined films considerably decreases in comparison with bulk 
liquids. This is manifested in the decrease of the diffusion constant and in the increase of the 
viscosity and molecular relaxation times. The "effective" viscosity can be 10 5 times larger 
and relaxation times can be 10 10 times slower than the bulk value. 

3. Boundary lubricant may exhibit two different responses in shear: a liquidlike, in which the 
liquid responds to the deformation by flow, and a solidlike characterized by observation of 
the development of yield stress in the confined fluid. In computer simulations with simple 
Lennard- Jones (LJ) liquids between atomically structured surfaces the existence of abrupt 
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liquid to solid transitions in films thinner than 6 molecular diameters was predicted, with 
the molecules becoming ordered both perpendicular and parallel to the surfaces. However, 
the horizontal order disappears for unstructured surfaces. In confined fluids both positional 
and orientational orderings are caused not only by mutual interactions of liquid molecules 
but also by the presence of two solid surfaces close together. 

4. For velocities of moving substrate less than some critical value the stick-slip motion may be 
observed. The stick-slip regime is characterized by intermittent stops (stick) and slips. 

Analytical models are widely used for explaining the experimental results. But these models 
are applicable only in special cases and, as a rule, they are only qualitative [ [H [H [10] . Computer 
simulations offer an additional tool for studying of the atomic-scale friction and wear. In particular, 
molecular dynamics (MD) simulations greatly help in the understanding of the nanotribological 

processes [ I El [12 US! E] • 

In this article the MD simulations of ultrathin water film confined between atomically flat 
diamond plates are described. The choice of water as a liquid for modelling is motivated by its 
importance, ubiquity and unique, in particular tribological properties [[3 [15]. An introducing of 
a water monolayer (only 0.25 nm thick) between atomically flat mica surfaces of the SFA causes 
the decrease of friction more than by the order of magnitude [ 16 . The effectiveness of a water 
monolayer to lower the friction force is attributed to the "hydrophilicity" of the mica surface 
(mica is "wetted" by water) and to the existence of a strongly repulsive short-range hydration 
force between such surfaces in aqueous solutions, which effectively removes the adhesion-controlled 
contribution to the friction force. Under the described conditions the friction force obeys the first 
Amontons' law which states that frictional force is proportional to the external load. 

The objective of the present work is to investigate the behaviour of the water film consisting 
of TIP4P molecules confined between totally rigid flat diamond plates and to check the model for 
consistency with the experiments. This in turn may give us estimates of the reliability of using the 
TIP4P model for water molecules, and the applicability of the approximation of absolutely rigid 
plates in the modelling of the tribological systems. 

2. Model 

The simulations are performed in a planar, Couette geometry that closely resembles the exper- 
imental systems [ IH [5] . A thin film of water molecules is confined between two solid walls with 
periodic boundary conditions in the plane of the walls. Each wall consists of 1152 atoms forming 
two (001) crystallographic planes of diamond lattice. Although the experimental study of friction 
of the diamond plates separated by the ultrathin water film has not been carried out yet but the 
diamond is transparent for the light and therefore can be used in the SFA. The walls are consid- 
ered as absolutely rigid and the elasticity of plates is not included into the model. But taking into 
account that diamond is one of the hardest materials we decided to verify such approximation. The 
water films of one and two molecular diameters thickness have been investigated. At the beginning 
of simulations one water layer consists of 196 molecules. The maximum value of particles involved 
in the simulations is 2696. 

Figure [I] shows the initial configurations of the studied system for both values of the film's 
thickness. The initial state of the film and of the crystalline plate (view is from the top in the 
negative direction of the axis z) is shown on figure [2] 

The empiric values of the covalent radii of atoms are used for visualization. For carbon, oxygen 
and hydrogen they are 77, 73, and 37 pm (1 pm = 10 -12 m), respectively, the lattice constant 
of diamond is 356.68 pm [[T7], At the beginning of simulations water molecules are positioned at 
the sites of a simple cubic lattice with the lattice constant corresponding to the density value of 
1048 kg m -3 . Molecular orientation is randomly assigned, with each linear and angular velocities 
having a fixed magnitude defined by the temperature and a randomly chosen directions. The water 
film and the plates are incommensurate because diamond lattice constant is not equal to the initial 
distance between water molecules. The initial distance (the gap) between the walls is 0.31 nm 
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Figure 1. The initial configuration of the studied system with one (left) and two (right) layers 
of water molecules. Carbon, oxygen and hydrogen atoms are displayed as blue, red and green 
balls, respectively. 
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Figure 2. The arrangement of water molecules (left) and the form of diamond plate (right) in 
the xy plane at the beginning of a simulation. 




for one layer of water molecules and 0.62 nm for two layers. The horizontal dimensions of the 
computational cell are equal for x and y directions and their value is 42.81 A. 

For imitation of the experiments the modelling is performed for constant temperature and load 
applied to the plates. The constraint method is used for maintaining the constant temperature 
value which is 298 K. The external load is modeled by applying the constant force L to each atom 
of the walls along the z axis. The constant-force algorithm is used [E]. This means that shear is 
imitated by applying constant horizontal shearing or driving force F$ to each atom of the upper 
plate along the x axis which corresponds to [010] direction. 

The TIP4P model is used for water molecules [ 18 . It is based on four interaction sites located 
in planar configuration, two of which - labeled M and O - are associated with the oxygen nucleus, 
and two - labeled H - with protons. The interaction energy between two molecules i and j consists 
of a double sum over the interaction sites of both molecules; the terms in the sum, indexed by k 
and /, allow for Coulomb interactions between the electric charges assigned to the sites, as well as 
contributions of LJ type, 




where r&/ = — r/ is the distance between sites k and /, rki = The corresponding force is 
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The charges appearing in the potential function are = 0.52e, qo = 0, <?m — — 2(/h, where e — 1.6 x 
10 — 19 C. The parameters in the LJ part of the potential, which acts only between O sites, are Aqq = 
A = 600 x 10 3 (kcal/mole)A 12 , C 00 = C = 610 (kcal/mole)A 6 , and 1 kcal/mole=4184 J/mole. 

The dimensionless units are used in the simulations. The length unit is a — (A/C) 1/6 = 3.154A, 
the energy unit is e — A/ (4<j) 12 = 0.155 kcal/mole, the mass unit is the mass of a water molecule, 
m = 2.987 x 10 _2T g and the time unit is to = \]m(j 2 js — 5.253 x 10 -12 s. Also in dimensionless 
units #h = 1- The time step used is At = 0.0005; in real units this equals to 2.627 x 10 -15 s. 

For the water-diamond interactions only the interactions of O site of the water molecule with 
carbon atoms are assumed. The interaction potential is of LJ type 



Uij= I A£C °[\^ ) \n 3 J j> '*J-'c ^ (3) 

[ 0, rij > r c 

where r c = 7.5 A(or 2.38 in dimensionless units) is the interaction cutoff distance, and sco — 2, 
ceo = 0.86. The values of latter parameters are chosen to fit the experimental fact that diamond 
is hydrophilic [ [19] . They correspond to the situation when the attractive forces between water 
molecules and carbon atoms have a twice greater magnitude than intermolecular water forces. The 
interaction between atoms of different plates is not introduced into the model. 

The classical equations of motion are used in simulations. As a consequence of the rigidity of 
diamond plates the motion of their centers of masses is concerned. The equations of motion for the 
upper plate are 

MX = F x + F s iV p , 

MY = F y , (4) 
MZ = F Z + LN P , 

where X, Y, Z are the coordinates of the center of mass of the upper plate, N p is the number of 
atoms of the wall, M = N p mc is the plate's mass, mc = 0.67 is the carbon atomic mass, F x , F y , F z 
are the components of the total force acting on the wall from the water. They are defined in a 
standard manner as the sum of the appropriate coordinate derivatives of the potential (3) with 
the negative sign for all molecules. Equations of motion for the lower plate are similar to (4). The 
difference is in that there is no component responsible for shear as there is in the first equation in 

@. 

The equations of the translational motion (without a component responsible for the constancy 
of temperature) for the center of mass of the ith molecule are as follow: 

ri = ^F ij + ^f ife , (5) 

j k 

where ^ij is the force acting from all carbon atoms that are not further that r c from the ith 
atom and it is defined in the way similar to the force for the surfaces; J2k ^ ^ s ^ ne f° rce acting on 
the ith molecule from other water molecules inside the cutoff distance, it is defined by the equation 
The rotational motion of water molecules is taken into account in the model and the equations 
of the rotational motion are expressed in terms of Hamilton's quaternions [HE]. 



3. Measurements, results and discussion 

Molecular dynamics simulations are performed for values of the load L per atom in the range 
from 2 (which is equal to 6.838 pN in real units) to 50 (0.171 nN). Corresponding pressure acting 
on the plates varies from 0.43 to 10.751 GPa. The horizontal shearing force F$ acting on each atom 
of the upper plate has values from 0.5 (or 1.71 pN) to 10 (34.3 pN) per atom and the resulting shear 
force acting on the plate changes from 1.97 to 394 nN. The maximum duration of the simulation 
is 52000 time steps or 136.6 ps. At each simulation during the first 2000 time steps the system is 
equilibrated and after that the measurements are carried out. There are two groups of simulations. 
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The first group is characterized by the absence of shearing force and the measurement of the 
diffusion constant is performed. In the second group the shear of the upper plate takes place and 
the friction force is measured. The latter is presented by the first component in Q. 

The diffusion constant D is calculated in two ways. In the first case the Einstein expression is 
used [HU 



D = lim 



i m 



fW-r^O)] 2 ), 



(6) 



where N m is the number of molecules, t is the time of measurements and Tj is the radius vector of 
the jth molecule center of mass. In ([6| angled brackets denote the average over a sufficiently large 
number of independent samples. 

The second way for calculating the diffusion constant is using the alternative Green-Kubo 
expression for D based on the integrated velocity autocorrelation function ip(t) [ [18| 120] 



D 



(7) 



(8) 



where Vj is the velocity of the jth molecule. 

Figure [3] shows calculated time dependencies of the velocity autocorrelation function of water 
molecules for different loads. It is apparent that with the increase of load the number and frequency 
of oscillations increases. This indicates that molecule velocities become more correlated under 
higher loads. 
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Figure 3. Time dependencies of the velocity autocorrelation function ^ for different loads and 
films with thickness of one (left) and two (right) molecular diameters. 

A gradual decrease of the diffusion constant with the external load increase shown on figure [4] 
and figure [5] indicates the transition of the water film to a solidlike state. The feature of this state is 
the ordering of molecules which can be vertical (in the form of quasidiscrete layers) and horizontal 
(or in-plane). Analytically the former manifests itself in the density oscillations with distance 
normal to the boundaries and the measure of the latter is the Debye- Waller factor. In this work 
mentioned quantities have not been calculated but ordering is observed visually. Figure [6] shows 
typical vertical ordering of water molecules for two values of load. In-plain ordering for one-layer film 
is not observed for all used values of load and driving force. However, for two molecular diameters 
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Figure 4. Time dependencies of the diffusion constant for different loads and films with thick- 
ness of one (left) and two (right) molecular diameters. Values are calculated with the Einstein 
expression ([61. 




Figure 5. The diffusion constant calculated with the Green-Kubo formula |7| as function of load 
for films with thickness of one and two molecular diameters. 



thick film in-plane ordering takes place especially for high loads. Typical molecule configurations 
are shown on figure [7[ Here one can see disordered and several types of ordered configurations of 
molecules. The presence of shear, as a rule, promotes the occurrence of the ordered structure for 
lower loads than without shearing. Also shear results in more rapid ordering. Thus, for two-layer 
thick water film there is a shear ordering that is observed experimentally [ H] . It should be noted 
that for low driving forces and intermediate loads molecular configuration is not completely defined 
by parameters of these forces. For example, structures shown on figure [7] for F$ = 0.5, L = 25 and 
Fs = 0, L = 30 are not absolutely reproducible and in other simulations for these values disordered 
configurations are obtained. In contrast to such behaviour for higher loads ordered state is always 
obtained as shown on the bottom right part of figure [7| 

The absence of ordering in the one-layer film we can explain after realizing the model peculiar- 
ities. The behaviour of one-layer film should be similar to the one of a simple LJ liquid confined 
between unstructured (i.e. mathematically smooth) surface. This follows from the facts that the 
roughness of surfaces is negligible and in such thin water film in the z direction its molecules can in- 
teract only with surface atoms. Surface atoms in turn interact only with O sites of water molecules 
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Figure 6. Layering of water molecules for values of load 8 (left) and 14 (right) and zero shear 
force. 

and only through the LJ potential. Therefore LJ contribution dominates in the interaction of water 
molecules in one-layer film. For the thicker film water molecules in adjacent layers can interact with 
each other and Coulomb interaction can cause the occurrence of ordered structure. It should be 
noted that observed ordering is not the true thermodynamic phase transition. And a liquidlike and 
a solidlike states are not really the same as the bulk liquid or solid phases. It would be more correct 
to refer to them perhaps as static and dynamic "epistates", since they arise only in boundary films, 
whose properties are determined both by the confinement and the epitaxial interactions between 
the film and surface atoms [0]. 

Figure [8] shows time dependencies of friction force obtained for different values of shearing force. 
For one-layer film such dependency is irregular with considerable fluctuations. For the thicker film 
periodic spikes are observed for higher values of shearing force. Similar spikes are observed in 
experiments when the film is in a solidlike state and they correspond to the stick-slip friction [ 
HJE]. But the considerable difference in time values in simulations and experiments which is 10 
orders of magnitude prevents the definite comparison of the results of simulations with experiments 
(in experiments time values are of order of seconds- minutes). 

Time averaged friction force as function of load is shown on figure [9j Here two main regions 
can be marked out. The first one is linear for low L and it corresponds to the first Amontons' law. 
The second part is horizontal and is observed for higher loads. Such dependencies can be explained 
with the use of the "cobblestone model" [ |3j [5j [16] . Accordingly to this model friction force is 
defined by two contributions. The first one arises from the presence of the internal adhesion forces 
between liquid and surface molecules. The second contribution is caused by the externally applied 
load. Frictional force is defined by the following expression 

F = S C A + CL, (9) 

where constants S c and C are critical shear stress and the friction coefficient, A is the area of 
contact and L is the external load. Quantity S c depends on the adhesion interactions of the film 
with surfaces. The friction coefficient C is related to the atomic granularity of the surfaces and on 
the size, shape, and configuration of the liquid molecules in the gap. In general, the smoother the 
surfaces the smaller should be the value of C. 

Dependencies on figure [9] may be explained in the following way. As the surfaces in simulations 
are smooth and water molecules are small and have a simple shape we can assume that frictional 
force is defined only by the adhesional contribution. Therefore we can consider that C « and in 
([9| there is only the first term in the right part of the equation. With the increase of load the number 
and strength of the adhesion bonds of water molecules with surfaces at first rapidly increase. This 
causes the increase of S c and from ([9| follows the increase of the friction force with L shown on 
figure [9j After reaching some external load for the current shearing force some stationary values 
of the number and strength of adhesion bonds are established. With the increase of L the contact 
area does not increase, since the model does not take into account the deformation of surfaces. So, 
the first term in ^ (and therefore friction force) remains constant. The same behaviour exhibits 
the shear stress S = F/A that is widely used in experiments. Note that in experiments for simple 
spherical molecules similar dependencies of shear stress S on load L have been obtained [15]. But 
for real surfaces the constancy of S for high L is caused by the proportional increase of F and 
A because of the deformation of surfaces. In contrast in simulations the constancy of S is the 
consequence of the constancy of both F and A. 
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Figure 7. Configurations of molecules in xy plane for different loads and shear forces at the end 
of corresponding simulations with the two-layer water film. 



Thus, considered model gives results in many ways corresponding to simple spherical molecules 
and reproduces for water under low loads the experimentally observed first Amontontons' law. 
Nevertheless, water molecules cannot be referred to as simple spherical. So, for obtaining more 
realistic results the model should be improved. For example, the elasticity, roughness, and mutual 
interactions of surfaces should be taken into account. Such more complicated models are the subject 
of the further investigations. 
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Figure 8. Time dependencies of the friction force for films with thickness of one (left) and two 
(right) molecular diameters and different driving forces. The value of load is 50. Negative trend 
is the consequence of the fact that friction force acts in the direction reverse to sliding. 
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Figure 9. Time-averaged value of the friction force for films with thickness of one (left) and two 
(right) molecular diameters and different shearing forces as function of load. 



4. Conclusion 

The main result of simulations is the transition of the ultrathin water film to a solidlike state. It 
is manifested in the decreasing of the diffusion constant and in the ordering of molecules. One layer 
of molecules between flat rigid diamond plates behaves like simple LJ boundary liquid lubricant 
and there is no horizontal ordering in it. For two-layer thick film the layering and in-plane ordering 
are observed. The presence of a driving force causes more rapid forming of ordered structures and 
the shear ordering is observed. Such behaviour can be explained by the peculiarities of the model, in 
particular by the presence of the Coulomb interactions between water molecules. Time dependencies 
of friction force in general reflect a solidlike state of the film and the stick-slip motion is observed 
for the thicker film under high shear force. But significant difference of time scales in simulations 
and experiments prevents the definite comparison of corresponding results. The dependencies of 
time averaged values of friction force and shear stress are similar to experimentally obtained for 
spherical molecules. They can be explained with the "cobblestone" model in the approximation of 
the adhesion forces domination. 
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